package edu.gmu.atelier;

/**
 * Performs QR factorization on the passed matrix. First executes the
 * GramSchmidt (need to change to house holder reflections) algorithm
 * to produce an orthogonal basis that spans the 
 * same space as A, then normalizes to produce an orthonormal basis.
 * <p>
 * Because (note only a left inverse) Q^T Q = I, we can solve for R in
 * the equation A = QR, such that R = Q^T A.
 * 
 * @author James H. Pope
 */
public class QRGramSchmidtFactor implements QRFactor
{
    private Matrix a  = null;
    private Matrix q  = null;
    private Matrix r  = null;
    
    public QRGramSchmidtFactor( Matrix a )
    {
        this.a = a;
        
        /*
         * Logically, after (i believe) n iterations, triangular orthogonalization
         * A R1 R2 ... Rn = Q
         * 
         * GramSchmidt has 2 drawbacks
         * 1. Requires A to have independent column vectors
         * 2. Has some numerical instabilities (though we do mitigate)
         * Not sure about square or rectangular, though (1) may imply?
         * 
         * Go look at http://www.math.iit.edu/~fass/477_577_handouts.html
         * Chapter 4 seems to provide a good explanation of QR Factor.
         */
        GramSchmidt gs = new GramSchmidt(a);
        //System.out.println("GS: Q="+gs.getQ());
        this.q = gs.getQ();
        this.r = q.transpose().mult( a );
    }
    
//    public static Matrix decompose( Matrix a )
//    {
//        Matrix r = a.copy();
//        
//        /*
//         * Using House Holder reflectors
//         */
//        List<Vector> rv = r.getVectors();
//        int n = r.getCols();
//        int m = r.getRows();
//        
//        // Remember the u's, they become the matrix Q
//        List<Vector> us = new ArrayList<Vector>();
//        
//        // Determine the maximum number of iterations required
//        // Note that this allows m=n, m<n and m>n matrices
//        int t = Math.min( m, n );
//        //int t = 1;
//        //int t = Math.min( m-1, n );
//        Matrix i = r.identity();
//        for (int k = 0; k < t-1; k++)
//        {
//            Matrix minor = r.copy(k, t-k);
//            System.out.println("Minor="+minor);
//            
//            Vector x = minor.getVector(0);
//            
//            Element kth = x.getElement(0);
//            Element norm = x.norm2().sqrt(); // this is definitely positive
//            if( kth.isNegative() )
//            {
//                norm = norm.negate();
//            }
//            
//            Matrix fi = minor.identity();
//            //System.out.println("Fi="+fi);
//            Vector e1 = fi.getVector(0);
//            
//            /*
//             * NOTE: This is not complete for Complex Matricies, need to test
//             * If real w always equals 1, therefore will be 2 times the vvt
//             * If complex, w = x^T.v / v^T.x
//             * I believe this supports both?
//             */
//            Vector u = x.sub( e1.mult(norm) );
//            us.add(u);
//            System.out.println("u="+u);
//            Vector v = u.normalize();
//            Element one = kth.getOne();
//            Element vtx = v.dot(x);
//            Element xtv = x.dot(v);
//            Element w = xtv.mult( vtx.inverse() );
//            Element multiplier = one.add(w);
//            //System.out.println("Multiplier="+multiplier);
//            Matrix twovvt = fi.mult(v, v).mult(multiplier);
//            
//            Matrix f = fi.sub(twovvt);
//            System.out.println("F="+f);
//            Matrix q = i.copyElements(k, k, f);
//            System.out.println("Q="+q);
//            
//            // First part of Q is the I identity matrix
//            r = q.mult(r);
//            System.out.println("R="+r);
//        }
//        return r;
//    }
    
    public Matrix getA()
    {
        return this.a;
    }
    
    public Matrix getQ()
    {
        return q;
    }
    
    public Matrix getR()
    {
        return r;
    }
}
